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R. C. Strawn and C. Tung 


U.S. Army Aeroflightdynamics Directorate, AVSCOM 
NASA Ames Research Center 
Moffett Field, California 94035 USA 


SUMMARY 

Two different schemes are presented for including the effect of rotor wakes on the finite-difference prediction of rotor 
loads. The first formulation includes wake effects by means of a blade-surface inflow specification. This approach is 
sufficiently simple to permit coupling of a full-potential finite-difference rotor code to a comprehensive integral model for 
the rotor wake and blade motion. The coupling involves a transfer of appropriate loads and inflow data between the 
two computer codes. Results are compared with experimental data for two advancing rotor cases. The second rotor- 
wake modeling scheme in this paper is a “split potential” formulation for computing unsteady blade-vortex interactions. 
Discrete vortex fields are introduced into a three-dimensional, conservative, full-potential rotor code. Computer predictions 
are compared with two experimental blade- vortex interaction cases. 

LIST OF SYMBOLS 

doo = freestream speed of sound 

AR = aspect ratio, R/c 

e = rotor chord length 

Ci = lift coefficient 

C p = pressure coefficient 

C T — rotor thrust coefficient 

G = prescribed velocity potential 

h = nondimensional time step 

1 = identity matrix 

J = Jacobian of coordinate transform matrix 

Mr — hover tip Mach number 
n = superscript representing current time level 

r = distance along the rotor in the spanwise direction 

f = radial coordinate in Fig. (1) 

R -- radial distance to the rotor tip 

Roo = free stream subtraction term in Eq. (5) 

t — normalized time in the untransformed coordinate system 

U — contravariant velocity normal to the plane 

V = contravariant velocity normal to the plane 

W = contravariant velocity normal to the {, r; plane 

x, y, z — untransformed coordinate system 

X = chordwise distance along the rotor 

i/ = spanwise coordinate in Fig. (1) 

= vortex core distance to the airfoil centerline 
a = angle of attack, deg. 

F = jump in potential across the wake normalized by (aooc) 

r„ = vortex strength normalized by (Oooc) 

"I = specific heat ratio 

li = rotor advance ratio 

$ = total velocity potential 

<j> — perturbation velocity potential 

i /i = rotor azimuthal angle, deg. 
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p = fluid density normalized by free-stream values 

a = rotor solidity 

r =■ normalized time in the transformed coordinate system 

P P/J 

0 ■— coordinate angle in Fig. (1) 

£, = transformed coordinate system 


1. INTRODUCTION 

Accurate methods for flow prediction on helicopter rotors are an important part of the rotor design process. Integral 
methods with prescribed wakes are routinely used for these flow-field predictions. An example of this type of approach is 
found in Ref. 1. However, integral techniques, such as panel and lifting-surface methods, are unable to treat transonic flow 
nonlineararities. This problem is a fundamental limitation of this methodology. 

Alternatives to integral methods arc finite-difference schemes that compute three-dimensional transonic flows around 
rotor blades. These methods typically solve the Euler equations, full-potential equations, or the small disturbance potential 
equations. These finite-difference methods have the advantage over other methods because they are designed to compute 
the transonic flow nonlinearities that are associated with high-speed advancing rotors. The disadvantage of these methods 
is that computer speed and storage requirements limit the size of the computational domain. The complexities of the rotor 
wake system usually extend well outside of typical finite-difference grids for helicopter rotors. 

To surmount these problems, several hybrid approaches have been developed that couple integral methods for the 
rotor wakes to finite-difference methods for the transonic flow over the rotor blades. Examples of these methods are given 
in Refs. 2 through 7. These approaches have been quite successful in predicting the aerodynamic loading on realistic 
helicopter rotors. 

Reference 7 presents an unsteady, three-dimensional, finite-difference method for computing the flow-field over ad- 
vancing helicopter rotors. This method solves the full-potential equations in strong conservation form. This paper also 
demonstrated a hybrid approach for including a rotor wake system into a hover computation. Vortex elements located 
close to the blade were modeled as discrete vortices inside the finite-difference computation. Vortices in the far-field were 
included using an integral wake model. 

The present paper has two objectives, both of which build on the results presented in Ref. 7. The first objective is 
to couple the Full- Potential Rotor code (FPR) in Ref. 7 to the integral wake modeling of the CAMRAD code in Ref. 1. 
The coupling scheme will be similar to those presented in Refs. 2 through 6. Although the coupling scheme has been used 
before with other finite-difference codes, the objective is to demonstrate its use with the FPR code. 

The second objective in this paper deals with a much more detailed representation of rotor blade-vortex interactions. 
The idea is to extend the finite-difference vortex modeling in Ref. 7 to compute unsteady problems. A single vortex is 
placed directly into the computational domain and the blade-vortex interaction is examined in detail. At present, the 
added complexity of the unsteady vortex limits this analysis to small numbers of vortices and relatively simple geometries. 


2. FORMULATION OF THE FULL-POTENTIAL EQUATIONS 


The basic FPR code used in this paper is described in Ref. 7. It solves the unsteady full-potential equation in strong 
conservation form. The rotor code is based on the fixed-wing code of Bridgeman et al. [8). 

The full-potential equation is written in a generalized coordinate system as follows: 


with density given by: 


dr V' + dC J ' + dv'' J ' + J ' 


P = {l + [-29> r ~(U + f t )* t - (V + .,()*„ ~(W + ft}#,] } 

where 


(1) 


(2) 


£/ = & + Ai$ { + A 4 $ n + A 6 $, 

V = i)t + A4${ + A 2 #,, + As#, (3) 

kP = ft + A 5 #f + A e #, + A 3 #; 


and 
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M = & + {* + f* 

A* = vl + nl + r)l 

Aa = £ + $ + C z ^ 

Ad — txVz 3" £y*ly 3" £zVz 

As ~ 3- £j/Cy 3- 

Aft ~ Vx$z 3" Vy$y "3* T}z$z 

All velocities are normalized by Ooo> distances by the airfoil chord length, and time by the combination (c/a rjci ). Density is 
normalized by the free-stream value. 

Equation (1) is solved by using first-order backward differencing in time and second-order central differencing in space. 
The temporal density derivative is locally linearized about the old time levels in a manner that preserves the conservative 
form (see Ref. 8). 

The resulting equation is approximately factored into r), and f operators. 


• 2 * 2 

[i + hu n s ( - X [/ + /.V% - 

( ' 2 r ( 5 ) 

[/ + hW n S s - -^6 ( (Ma)%] (* n+l - $ n ) = ^ [h(pV) n + Sy(pV) n + S s (pW) n - Rn] + C 

where S(, S^, and 6$ represent central difference operators in space and the term C is given by Eq. (6). The bracketed 
term in Eq. (6) represents the temporal conservation correction to the algorithm. 


C =($ n - $ n -M + - 2$ n ~ 1 + $"- 2 ) + i-(p n - p"- 1 

L f)n 

+ hjr^-{U n - 1 6 ( + K"- 1 ^ 3- W' n " 1 5 f )($ n - S”- 1 )] 

A steady state ADI relaxation algorithm can be obtained from Eq. (5) by omitting the unsteady C term on the 
right-hand side of the equation. 

The streamwise flux terms use upwind density biasing in regions of supercritical flow to ensure stability of the algorithm. 
Details are given in Ref. 8. The quantity R&, represents a numerical truncation error term caused by incomplete metric 
cancellation. Formulation of this correction term is discussed in Ref. 8. 



3. GRID AND BOUNDARY CONDITIONS 

Rotor flow fields are computed by modifying the time-metric terms t/ t , and ft in Eqs. (2) and (3). The finite- 
difference grid is attached to the rotor blade which moves through still air. Each grid point has an appropriate rotational 
and translational coordinate velocity field given by £ t , r) t , and ft- This results in free-stream conditions that are given by 
$ = 0. A spanwise series of planar parallel O-grids have been chosen for the finite-difference mesh. This results in a very 
efficient use of grid points. Grid size for a typical rotor case is 80 chordwise points, 19 spanwise points, and 25 points 
normal to the surface. The radial location of the outer grid boundary typically varies from 7 to 10 chords from the rotor 
surface. 

A surface transpiration condition is implemented on the rotor surface in order to simulate angle of attack conditions. 
The transpiration condition can be written as ~ ($„ + St) tan a, where a is the desired angle of attack, s t represents the 
surface coordinate motion, and the s and n directions are tangential and normal to the airfoil surface respectively. Along 
the inner boundary plane, normal to the rotor, the spanwise contravariant velocity V, is set equal to tje 

For lifting cases, the shed vorticity is specified as a jump in potential, T, imposed across the coordinate plane £ = 0. 
This coordinate plane is approximately aligned with the shear layer from the trailing edge. An unsteady transport equation 
for this potential jump is imposed across the wake. It is derived by assuming that density is continuous across the cut, and 
can be written as: 


r r + <v>r„ + qv)T t = o (7) 

where ( V ) and (IF) are the averages of the contravariant velocities above and below the wake. 

At the outer grid boundary, a nonreflection boundary condition has been implemented similar to that used in Ref. 9. 
If we make the assumption that disturbances propagate only in the radial and spanwise directions at the outer boundary 
of the O-grid, then a nonreflection boundary condition can be written as: 


<t>t + ((Moo)* COS 0 -l)0f + [Moo)y4>y = 0 


( 8 ) 
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where the radial coordinate f and the angle 6 are defined in Fig. 1. The values (Moo)* and are the free-stream 

Mach number components relative to the blade-fixed grid. 

Equation (8) is used to explicitly update the solution on the outer boundary. It is applied after each time-step for 
the interior points. Implementation of Eq. (8) was necessary for unsteady lifting rotor cases so that the accumulation of 
disturbances at the outer boundary could be prevented. 

4. WAKE MODELING FOR FORWARD FLIGHT COMPUTATIONS 


Lifting rotors in forward flight have wake systems that are too complicated to model directly in the current finite- 
difference full-potential scheme. In addition, there is no easy way to couple the blade blade trim solution directly into a 
finite-difference model. A way to circumvent these difficulties is to couple the FPR code for the rotor aerodynamic loads to 
the CAMRAD integral method [1] for modeling the rotor wake and blade motion. The basic coupling scheme is the same 
as that described in Refs. 4 and 5. 

The CAMRAD code is a fast and efficient model for comprehensive rotor performance analysis. It normally models the 
aerodynamic loads through equivalent angles of attack for the rotor that result from the trimmed-rotor and wake system. 
It then uses table lookups of two-dimensional airfoil data to determine the blade lift coefficients. These table lookups do 
not model the flow-field unsteadiness or three-dimensionality that is present in actual advancing rotor cases. The coupling 
scheme used in this paper replaces these table lookups with lift coefficients that are determined from finite-difference 
calculations for the specified partial angles of attack. 

The matching of the CAMRAD and FPR codes is summarized in Fig. 2. The process begins by obtaining a trimmed 
solution from the CAMRAD code with the lift totally obtained from airfoil tables. The CAMRAD solution specifies 
equivalent angles of attack along the rotor. This rotor twist distribution is used to model the wake effects in the subsequent 
FPR finite-difference calculation. Unsteady lift values resulting from the FPR calculation are then applied to the next 
CAMRAD trim solution using the following algorithm. 

Ci(a) = Ci PPR (ctoid) + C llabu (a) - C ltMe {a M ) (9) 

The blade twist angles a and a 0 u are the angles of attack from the current and previous trim loops, respectively. The 
particular form of Eq. (9) is chosen to enable the CAMRAD computer program to dynamically update the lift coefficient 
solution during the calculation. The solution is converged when a = a 0 id and the finite-difference lift is fully consistent 
with the rotor inflow and motion. Convergence of the scheme is rapid and typically occurs in less than three iterative cycles 
around the loop in Fig. 2. 

5. RESULTS FOR ADVANCING ROTORS 

The CAMRAD and FPR code coupling was tested in two cases. The first case is a 1/7-scale model Cobra operational 
load survey (OLS) rotor that is described in Ref. 10. This stiff-bladed AH-1G model was pressure-instrumented. Data is 
currently available only on the upper surface of the rotor at r/R — 0.95. Tip Mach number, Mr, is 0.663, and advance 
ratio, n, is 0.298. The finite-difference computation used 720 equal time-steps in advancing from 0° to 180°. Total computer 
time for each global iteration of the FPR code was approximately 10 CPU minutes on a CRAY XMP-12 computer. 

Figure 3 shows the CAMRAD/FPR pressure predictions for this case at a series of rotor azimuthal angles. Experimental 
data on the upper blade surface are also shown. Predicted results in the first quadrant show good agreement with the 
exerimental results. At tf> = 90°, the prediction shows more expansion at the leading edge than is present in the data 
although the shock location is predicted quite well. 

The negative surface pressures are overpredicted in the second quadrant when they are compared to the data. This is 
consistent with the results of Tung and Chang [4], who solved a nonconservative form of the full-potential equation coupled 
to the CAMRAD wake model. It is difficult to determine the exact cause of the discrepancy in the second quadrant. This 
is particularly curious in light of the fact that the wake structure in the first quadrant is more complex than that in the 
second quadrant. 

The second computed test case was a three-bladed rotor tested at ONERA by Philippe and Chattot [11], The blade 
was pressure-instrumented at three radial locations, r/R — 0.85,0.90, and 0.95. Both upper- and lower-surface pressures 
were reported in the experiment. Blade tip Mach number Mj- is 0.629, and the advance ratio is relatively high at n = 0.387. 
Once again, 720 equal time-steps were used from j/j = 0° to $ = 180°. 

Figures 4, 5 and 6 show the comparison between prediction and experiment at r/R — 0.85,0.90 and 0.95 respectively. 
At each radial location, pressure results on the rotor are presented for a series of azimuthal angles. Agreement between 
the prediction and the data is good, particularly for the upper surface shock location, and for the lower surface expansion 
near the leading edge. Predicted pressures on the lower surface are slightly higher than those shown by the experimental 
data, particularly in the second quadrant. However, predicted pressures on the upper surface, are in good agreement with 
the data. Overall, the evolution of the pressures is well predicted by the combined CAMRAD/FPR method. 



«. UNSTEADY BLADE- VORTEX INTERACTION MODEL 
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Reference 7 describes a method by which discrete vortices could be introduced into the finite-difference computation. 
The approach is to consider the difference equation (5) in functional form as shown below. 


L(p n ,$ r, ,$ n+1 - $") = R(p",$ n ) + C(pV n - 1 ) $\$ n -\$ n - 2 ) 
p" = F($ n ) 


( 10 ) 


The first relation represents mass conservation and the second relation is the Bernoulli equation. Note that the L 
operator in Eq. (10) is linear with respect to 0 n+1 - 4>". This is a result of the original density linearization. 

The total potential $ is now rewritten as the sum of two parts. The first is a known potential, G, which is due to a 
specified system of vortex elements. The second component is an unknown potential, rf> % which is the perturbation of G 
that is due to the rotor and its immediate shed vorticity. If we now substitute the equality, $ = G + <j>, into Eq. (10), the 
following form results: 



(p\P n - x \ 



L(p n ,# n ,<T +1 - <f> n ) =R (p n ,<t> n + G n ) + C 

r 

r~' 

+ c 

G" 

G n-1 


U "- 2 J 


tG" -2 7 


p n =F(d>” + G n ) 


-L (p n ,$ n ,G n+l -G n ) 


( 11 ) 


The left-hand side of Eq. (11) is identical to the original algorithm. The right-hand side and the Bernoulli relation undergo 
only minor modifications. 


For the steady hover problem considered in Ref. 7, the right-hand functions in Eq. (11) contains only spatial 
gradient terms in G. They can be added to the original FPR code with relatively little effort. An unsteady, forward-flight 
computation requires that the temporal gradient terms involving G must also be added. Because of the added complexity 
in programming, some of the temporal gradient terms involving G have not been included at present. The unsteady vortex 
terms that have been dropped from the right-hand side of the final solution scheme are 


(p n ,p n ~ 1 \ 

G n 

Q n-1 

G n ~ 2 


- L(/>",$ n ,G n+1 - G n ) 




J 


( 12 ) 


However, all spatial and temporal gradient terms resulting from the vortex field have been included in the density and 
contravariant velocity evaluations. 

If the center of a vortex is located far from the rotor surface, say greater than one chord length, then it can be shown 
analytically that the neglected terms in Eq. (12) are small. For vortex passage distances that are less than one chord, it 
becomes more difficult to justify this assumption. In this paper, typical vortex passage distances are approximately half a 
chord length. 

Jones [12] has studied the unsteady airfoil-vortex interaction problem for two-dimensional cases. He has modeled an 
unsteady vortex in a number of different ways. Vortex strengths and passage distances are similar to the three-dimensional 
cases considered in the present paper. Jones’ results show that pressure and lift results from a “branch-cut” vortex model 
do not differ appreciably from those of a “split-potential” method which neglects the terms in Eq. (12). In the branch-cut 
model, the vortex is included in the flow as a discontinuity in total potential. The branch-cut method implicitly includes 
all of the temporal gradient terms that are caused by the vortex, but in a way that is restricted to an H-mesh. References 
3 and 12 give a further description of the branch-cut vortex model. 

Sankar and Malone [13] also present a series of two-dimensional, blade-vortex calculations. Their model uses a split- 
potential formulation similar to the one described above. They neglect, however, all of the temporal and many of the 
spatial derivative terms that result from the vortex field. Vortex strengths and passage distances are similar to the cases 
considered in this paper. Their airfoil lift results compare favorably with the two-dimensional transonic small-disturbance 
results of McCroskey [14], the Euler code results of Sankar and Tang [15], and the full-potential results of Jones [12]. 

The evidence presented above, coupled with our own comparisons to both experimental data and other computer 
codes, indicates that the neglected terms in Eq. (12) are of little consequence for the vortex strengths and passage 
distances considered in the present paper. In future work we plan to add the neglected terms to the FPR code to verify 
this assumption. 

Strictly speaking, the vortex velocity field in these calculations must be irrotational since it is derived from the gradient 
of the potential G. In practice however, we include a small rotational core in our vortex velocity fields in order to prevent 
the velocities from becoming infinite at the vortex centers. 
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7. RESULTS FROM UNSTEADY BLADE- VORTEX INTERACTION CASES 


To test the vortex model described above, the FPR code was set up to compute the unsteady blade-vortex interaction 
experiment of Caradonna et al. [16]. In this experiment, a two-bladed NACA 0012 model rotor was was tested in a wind 
tunnel. A NACA 0015 wing was located upstream in the tunnel. It was used to produce a streamwise line-vortex that 
passed beneath the helicopter rotor. Figure 7 is taken from Ref. 17 and illustrates the experimental design. Pressure data 
was measured at a single spanwise location on the blade (r/R = 0.893). Blade aspect ratio was equal to 7. Srinivasan et 
al. [17] have previously presented computational results for this experiment using two-dimensional thin-layer Navier Stokes 
and transonic small-disturbance models. 

The first test case is entirely subcritical, with a tip Mach number of 0.6, and an advance ratio of 0.2 . At the 0° and 
180° azimuthal rotor angles, the line-vortex is located directly under the quarter-chord axis of the rotor blade at a distance 
of 0.4 chords. The vortex core radius for the calculation was set equal to the experimentally estimated value of 0.15 chords. 
Predicted pressure results showed little sensitivity to vortex core size. 

The blade- vortex interaction calculations were started with a quasi-steady solution at ij> = 30°. Each time-step covered 

0.25° of rotor azimuthal angle until the rotor reached 4> — 175°. At this point, the rotor is close to the vortex and the 
time-step was cut to 0.1°. Computer time for the calculation is approximately 0.83 seconds per time-step on a CRAY 
XMP-12 computer system. Total CPU time for these calculations was approximately 10 minutes. 

There is some question as to the accuracy of the vortex strengths reported for this experiment. Srinivasan et al. 
[17] have compiled a large amount of experimental data for vortex velocity distributions from NACA 0015 wings. Based 
on this data, their calculation used a nondimensional vortex strength T„ = 0.133, which was approximately 80% of the 
value reported by Caradonna et al. [16]. The calculations in this paper use the vortex strengths that were determined by 
Srinivasan et al. [17]. 

Figure 8 shows the results from this calculation. Pressure coefficients are presented at the r/R = 0.893 spanwise station 
for a series of rotor azimuthal angles. The predictions show excellent agreement with the experimental data, especially 
in the crossover region from negative-to-positive lift. The main discrepancy between the prediction and the data is the 
presence of a slight pressure wave in the predicted lower surface pressure. This pressure wave is located directly above the 
vortex as the rotor passes over it. It does not appear to be present in the experimental data and may be the result of 
neglecting the temporal derivative terms in Eq. (12). 

Computations have also been performed for a supercritical case, Mr = 0.8. The data in this case involves interactions 
where the shock persists up to and throughout the interaction. To date, all computational attempts at this case have 
been unsuccessful because the computed shock collapses prematurely in the second quadrant (by approximately 5°). The 
reasons for this are not yet known. 

To demonstrate that the FPR code can qualitatively predict a transonic blade-vortex interaction if the shock persists, 
the tip Mach number for the computer simulation was raised to Mr = 0.82. The results are shown in Fig. 9. The higher 
tip Mach number ensured that the shock structure was in better agreement with the data at the start of the blade-vortex 
encounter (xj> — 170°). With this adjustment, it is clear that the computer prediction demonstrates the major features of 
the transonic blade-vortex interaction. The shock on the upper surface disappears rapidly, while the lower surface shock 
moves forward more slowly toward the the leading edge of the rotor. The unsteady shock motion on the vortex side of the 
rotor is the dominant feature in the flow. It could have major effects on rotor acoustics and on the blade higher harmonic 
loading. 


8. CONCLUSIONS 


1. The FPR code has been coupled to the CAMRAD wake model. The result is a hybrid integral-differential rotor flow 

scheme. Successful prediction of two advancing rotor cases has been demonstrated. 

2. A “split-potential” formulation has been developed for computing three-dimensional, unsteady blade-vortex interactions. 

Discrete vortices are implicitly coupled to a conservative full-potential model for the rotor aerodynamics. 

S. The split-potential FPR code has been used to obtain good results for a subcritical blade-vortex interaction case. 
Qualitatively realistic results have been shown for a blade-vortex interaction in a supercritical case. Results thus far 
have been limited to simplified geometries and a single vortex. More complex blade-vortex interactions will be studied 
in future work. 
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Fig. 1 Coordinate system for a nonreflection boundary condition. 



DONE 


Fig. 2 A hybrid integral-differential rotor flow scheme. 
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UPPER SURFACE PREDICTION 

LOWER SURFACE PREDICTION 
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o LOWER SURFACE EXPERIMENTAL DATA 










Fig. 5 Surface pressure results for the ONERA three-bladed rotor, Mt = 0.629, n = 0.387, AR = 6.968, r/R = 0.90, 
Ct/o = 0.0665. 




7-13 



Fig. 7 Rotor-vortex interaction experiment in Ref. 16. 

LOWER SURFACE PREDICTION 

UPPER SURFACE PREDICTION 

• LOWER SURFACE EXPERIMENTAL DATA 
o UPPER SURFACE EXPERIMENTAL DATA 






Fig. 8 Surface pressure results for a three-dimensional blade-vortex interaction, = 0.133, Z„/c = 0.4, Mt — 0.6, 
A = 0.2, AR = 7.0, r/R = 0.893, a = 0, untwisted, untapered, NACA 0012 blade. 



7-14 


LOWER SURFACE PREDICTION (M T = 0.82) 

UPPER SURFACE PREDICTION (M T = 0.82) 

• LOWER SURFACE EXPERIMENTAL DATA (M T = 0.8) 
o UPPER SURFACE EXPERIMENTAL DATA (Mj = 0.8) 




Fig. 9 Surface pressure results for a three-dimensional blade-vortex interaction, F„ = 0.177, Z v /c = 0.4, p = 0.2, AR = 7.0, 
r/R = 0.893, a = 0, untwisted, untapered, NACA 0012 blade. Data has Mt = 0.8, prediction is for Mr = 0.82. 
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